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A continuum-level, dual internal state variable, thermodynamically based, work potential model, Schapery 
Theory, is used capture the effects of two matrix damage mechanisms in a fiber-reinforced laminated com- 
posite: microdamage and transverse cracking. Matrix microdamage accrues primarily in the form of shear 
microcracks between the fibers of the composite. Whereas, larger transverse matrix cracks typically span the 
thickness of a lamina and run parallel to the fibers. Schapery Theory uses the energy potential required to 
advance structural changes, associated with the damage mechanisms, to govern damage growth through a set 
of internal state variables. These state variables are used to quantify the stiffness degradation resulting from 
damage growth. The transverse and shear stiffness’ of the lamina are related to the internal state variables 
through a set of measurable damage functions. Additionally, the damage variables for a given strain state can 
be calculated from a set of evolution equations. These evolution equations and damage functions are imple- 
mented into the finite element method and used to govern the constitutive response of the material points in 
the model. Additionally, an axial failure criterion is included in the model. The response of a center-notched, 
buffer strip-stiffened panel subjected to uniaxial tension is investigated and results are compared to experi- 
ment. 


I. Introduction 

D amage and failure analysis of composite materials continues to be prevalent and challenging research topic. The 
intricacies of damage and failure in composite materials are attributed to numerous sources. The properties of 
constituent materials, the architecture of the reinforcement material, the interfacial behavior, and the applied loading 
configuration all influence the type and evolution of damage mechanisms in composite materials. Furthermore, many 
globally observed damage mechanisms are a result of the interaction of the damage mechanisms in the different phases 
and across various length scales. It is the interaction of damage mechanisms, due to the heterogeneity of the material, 
that results in composites behaving more like a structure and is what separates them from monolithic materials. In 
order to capture the non-linear behavior of composites correctly, the fundamental mechanisms must be modeled using 
the underlying physics behind those mechanisms. Approaching composite modeling in this manner can lead to very 
accurate, predictive models; which in turn, can lead to great weight and cost savings in design of composite structures. 

For unidirectional, fiber-reinforced, epoxy laminates, numerous damage mechanisms are observed.^ Damage first 
begins in the matrix of the composite at a scale on the order of the fibers. At low stress levels, matrix microdamage 
arises in the form of micro-void growth, shear microcracks and shear bands between fibers.^’ ^ This damage mechanism 
is highly distributed and is the primary cause for the globally observed non-linear behavior up to the onset of more 
severe damage mechanisms. As the load increases, transverse cracks begin to nucleate. These cracks typically run 
parallel to the fibers and span the thickness of a lamina in the fiber-reinforced laminate (FRL).^“^^ This mechanism 
occurs on a larger length scale than microdamage. Additionally, the effects of this mechanism on the overall stiffness 
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of the composite are more pronounced. In angle-ply and unidirectional coupons, the initiation of transverse cracking 
quickly leads to catastrophic failure of the specimen. Although these damage mechanisms are exclusive to the matrix 
of the FRL, they are heavily influenced by the flber properties and architectures. Placement of the flbers creates local 
stress concentrations that nucleate and drive cracks within the materials. Also, the stiff flbers are responsible for 
arresting matrix crack growth. Thus, the behavior of neat resin is different than the in-situ matrix 

Axial (flber-direction) damage of a ply in a laminate is attributed to damage in the flbers; however, the matrix 
and flber/matrix interface also play a role in axial damage. Carbon flbers, used in common FRL applications, usually 
behave linearly until the flbers rupture suddenly. However in tension, axial failure of composite plies within a laminate 
is often a combination of flber breakage and the flber being pulled out of the surrounding matrix. Fiber-pullout 
can occur as macroscopic cracks progress within a laminate. It can result when there is a through-thickness matrix 
crack, but only some, or none, of the flbers in that region have broken. Fiber-pullout can also arise when flber 
failure develops away from the primary fracture plane. In these scenarios the flbers must debond and “pull out” of 
the surrounding matrix in order for the crack to open further. Thus, flber-pullout is heavily dependent on the matrix 
properties and interfacial properties. The onset of axial damage typically indicates that catastrophic failure of the 
laminate is forthcoming; however the process can be slowed due to flber-bridging between the two crack faces. 

In compression, the globally observed axial damage mechanism is flber micro-buckling. Small misalignments 
of the flbers lead to flber rotation as the structure is loaded. These rotations induce matrix softening in shear; which, 
in turn, leads to further flber rotations. Eventually, this leads to a runaway instability and a flber kink band forms. 

Often, in-plane lamina damage is accompanied by delamination between the layers. Interlamina damage occurs 
when adjacent layers displace relative to each other and the bond between those layers is weakened. Delamination 
can result from edge effects, or it may develop when a matrix crack, within a ply, reaches the interface of the adjacent 
lamina. It has been shown experimentally that delamination can influence the globally observed failure modes of 
open-hole tension, composite specimens. At the same time, in compression of open hole specimens, it has been 
shown that flber kink banding is dominant and leads to delamination formation.^^"^^ 

There exist numerous approaches to modeling damage mechanisms in a composite. Some methods for capturing 
in-plane damage mechanisms at the continuum-level include using empirical failure criteria, continuum damage me- 
chanics, and/or fracture mechanics and discontinuous methods. In-plane damage mechanisms can also be 
captured with micromechanics models including repeating unit cell (RUC) and representative volume element (RVE) 
techniques, as well as other methods.^’ Moreover, micromechanics models can be tied to the continuum scale 
using multiscale methods in which the RUC is homogenized to give the effective response of the continuum.'^^"'^^ 
Commonly used techniques for introducing interlamina damage into a composite simulation include the virtual crack 
closure technique (VCCT), cohesive zone model and discrete cohesive zone model (DCZM).^^“^^ Ad- 

ditionally, multiscale methods can also be used to simulate interlamina effects. 

The method presented in this work uses three different models to encapsulate the effects of various damage mech- 
anisms in an FRL. Matrix microdamage and transverse cracking are modeled using a dual internal state variable (ISV) 
formulation of the thermodynamically based work potential model developed by Schapery, and presented in Ref. 24. 
Single and dual ISV formulations of Schapery Theory (ST) have been used to model microdamage alone and both 
microdamage and transverse cracking in analytical and numerical frameworks.^’ Lamina-level ST represents 
each layer as a homogenized continuum and uses the amount of energy necessary to advance structural changes within 
the material, associated with the damage mechanisms, to govern damage growth in a thermodynamically consistent 
manner. Each damage mechanism is accounted for through separate IS Vs which are related to the degraded stiffness’ 
through a set of measurable damage functions. In this paper, ST is implemented into the flnite element method (FEM) 
and is used to control the constitutive behavior of the elements. Additionally, a simple maximum strain axial failure 
criterion in used for axial failure. Upon satisfaction of the criterion, the stiffness of the corresponding material point 
is signiflcantly degraded. A more physically-based model for axial failure is desired and will be included in future 
models. Einally, delamination is included in the model via cohesive elements. The details of the flnite element model, 
including the ST implementation, and axial failure are presented in Section IV. 

Center-notched, IM7/8552 panels were tested in uniaxial tension at the NASA Langley Research Center.^^ Two, 
one inch sections, spanning the length of the panels, were reinforced with additional 0° plies. These sections, referred 
to as buffer strips, were included to increase the ultimate load of the panels, without introducing too much additional 
weight, and to change the failure mode of the panels. The experimental program is outlined in Section III. Section IV 
provides details of an EEM model of the buffer strip panel. The global load versus displacement and local strain gage 
data from the experiment are compared to the results of the EEM model in Section V. Additionally, analysis of the 
damage and failure mechanisms observed in the EEM model is provided in Section V. 
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II. Thermodynamically Based Work Potential Model 


A thermodynamically based framework for modeling elastic solids undergoing structural changes was developed 
by Schapery.^"^ In this theory the total work potential Wt is partitioned into a recoverable, elastic portion W, and a 
dissipated, or inelastic, part Ws- 

Wt = W + Ws ( 1 ) 

As the material is loaded, a portion of the applied work potential facilitates structural changes in the material. The evo- 
lution of these structural changes affects the elastic properties of the material. However, a fraction of the total applied 
work potential is recovered when the material is unloaded. The magnitude of energy recovered is contingent upon the 
current, degraded elastic properties. It is assumed upon subsequent reloading, the material behaves linearly, exhibiting 
the elastic properties observed during unloading, until it reaches the previous maximum strain state. After this state 
is achieved, structural changes resume further degrading the elastic moduli of the material. This secant assumption 
is not mandatory and may be relaxed depending on the material, but it is acceptable for fiber reinforced composites 
exhibiting low levels of plastic deformation and viscoelastic behavior.^ This process is displayed in Figure 1 in the 
form of a typical stress- strain curve for a material under uniaxial tension. The area under the linear unloading curve 
is the recoverable, elastic strain energy density W and the remaining, shaded area under the curve is the dissipated 
energy per unit initial volume Ws . 



Figure 1. Typical stress-strain curve displaying the elastic {W) and irrecoverable (tvs') potentials. 


Both W and Ws are functions of a set of internal state variables (IS Vs), Sm, = 1, 2, M). These IS Vs ac- 
count for any inelastic structural changes in the material. Differentiating Ws with respect to any ISV, Sm yields the 
thermodynamic driving force, /^, available for producing structural changes associated with the ISV. 


^ dWs 
dSm 


( 2 ) 


It is shown in Ref. 24 that the total work potential is stationary with respect to each ISV. 


8Wt 


dSrr 


= 0 


( 3 ) 


Additionally, Rice^^ has shown that according to the second law of thermodynamics: 


fmSm > 0 


( 4 ) 


Equations (2), (3), and (4) form the foundation of a thermodynamically based work potential model for nonlinear 
structural changes in a material. If one of the IS Vs is chosen to be crack length. Equation 3 is analgous to Griffith’s 
criterion commonly used in fracture mechanics. 


II. A. Dual-ISV formulation of ST to model progressive microdamage and transverse cracking 

The generality of the evolution equations, presented in the previous section, allows for the response of a material 
undergoing any number and type of structural changes to be captured. This is especially useful for modeling progres- 
sive damage in composites because multiple damage mechanisms arise during a typical loading history. For instance 
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in the matrix phase alone, microdamage accrues until its effects are superseded by the growth of larger transverse 
cracks. Microdamage is considered the advancement of microcracks, voids, fissures, shear bands, and other flaws 
that are present in the matrix of a composite.^’ The size of these fiaws is typically on the order of that of the 
fiber or smaller. Transverse cracks nucleate from preexisting fiaws within the matrix but grow parallel to the fibers 
and span the thickness of the lamina.^’ Often, the growth of individual transverse cracks is extremely 
rapid; however, the effects of transverse cracking on the stiffness of a composite laminate can be progressive if mul- 
tiple cracks form over an extended period of time. Eventually, transverse cracking is succeeded by more catastrophic 
damage mechanisms like interlamina delamination, fiber breakage, pull-out and bridging associated with macroscopic 
laminate fracture . 1 1 ’ 47, 6 1 

The average behavior of numerous coupon level laminates exhibiting progressive microdamage and transverse 
cracking was successfully modeled in Refs. 2 and 54 using ST. A similar formulation will be used in this work, in 
conjunction with the FEM, to model more complicated laminates. To incorporate microdamage and transverse crack- 
ing into ST, it is assumed that the inelastic work potential Ws is a function of two IS Vs: S which accounts for energy 
dissipation via microdamage, and Sc which captures available energy to form and advance transverse cracks. It can 
be assumed that these IS Vs are equal to the energy density required to advance the associated mechanism. Although 
interaction between the mechanisms occurs through coupling in the evolution laws and adjustment of the stress and 
strain fields with increasing damage, it is assumed that these IS Vs are independent; in that, microdamage evolution 
does not directly influence the growth of transverse cracks and vice versa. Thus, the total amount of dissipated energy 
is a sum of the available energy used to progress each individual mechanism: Ws = S' + 5'c- The total energy is then 
given by 

Wt = W^S^Sc (5) 


Applying stationarity of potential energy (Equation (3)) to the definition of the total work potential yields two 
evolution equations. 


dW 

dW 


-1 

( 6 ) 

-1 

(7) 


Additionally, invoking Equation (4) results in statements of the inadmissability of healing for both microdamage and 
transverse cracking. 


S > 0 (8) 

Sc > 0 (9) 


Equations ( 8 ) and (9) dictate that the amount of energy used to progress microdamage and transverse cracking can 
never decrease; therefore, that energy has been dissipated into creating structural changes and cannot be recovered. 
The combination of Equations ( 6 )- (9) represent the evolution equations for microdamage and transverse cracking in 
the matrix of the composite. 


II.B. Application of dual-ISV ST to fiber reinforced laminae 

II.B.l. 2-D plane stress constitutive law 

To tailor the evolution equations presented in Section II. A to a unidirectional, fiber reinforced lamina, the stress-strain 
relationships can be defined, under plane stress assumptions, in the principal material coordinates of the lamina as 

= QllCii -1- Qi2C22 

(^22 = Q \ 2 ^\\ + Q 22^22 ( 10 ) 

7"12 = Q66712 

where Cij are the stress components, are the strain components, the 1 -direction is aligned with the fibers in the 
lamina, the ’ 2 -direction is perpendicular to the fibers, 712 is the engineering shear strain and Qij are the components 
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of the stiffness matrix given by 


Qii 

Q22 


Ell 

1 -U12U21 
E22 

1 - Z/12Z/21 


Qi 2 

^66 

Z ^21 


— ^ 12^22 
= Gi 2 
_ 1^12^22 
Ell 


( 11 ) 


where En is the axial elastic modulus, E22 is the transverse elastic modulus, z/12 is the Poisson’s Ratio, z/21 is the 
transverse Poisson’s Ratio and G12 is the elastic shear modulus. After assuming that the quantity 1^121^21 « 1 , 
Equations ( 11 ) become. 


Qii = Ell 
Q22 = E22 
Qi2 = ^ 12^22 
Qee = Gi 2 


( 12 ) 


1LB.2. Determining the damage state 

It is necessary to define the manner in which the moduli degrade as functions of the IS Vs. Since the damage mecha- 
nisms associated with S and Sc accrue only in the matrix of the composite, it is assumed that the moduli affected by 
this damage are limited to E22 and G12. These moduli can be written as functions of S and S'c- 


E22 = E22oes{S)ec(Sc) 

( 13 ) 

Gi 2 = Guogs{S)gc{Sc) 

( 14 ) 


where E220 and G120 are the undamaged transverse and shear elastic moduli, 6 ^( 5 ') and gs{S) are factors relating 
the transverse and shear moduli to the microdamage S; edSc) and gdSc) are functions relating the moduli to the 
transverse cracking ISV Sc- Sicking provided a procedure for determining the damage functions experimentally.^ The 
experimental damage curves can then be fit with polynomials (such that values at S' = 0 and Sc = 0 , are E22 = E220 
and Gi 2 = G120) and used in Equations ( 13 ) and ( 14 ). 

The elastic strain energy density W 

= y cTijdeij ( 15 ) 

can be written using the plane stress constitutive relationships. 

w = -{Eii€ii + E22^22 + G12712) + ^12^11^22 ( 16 ) 

Employing Equations ( 6 ) and ( 7 ) with ( 13 ), ( 14 ), and ( 16 ), while assuming the quantity Q12 = ai2Q22 (the effects 
of evolving microdamage on the Poisson’s ratio have not be studied in detail) is constant (independent of S and Sc), 
yields two ordinary differential equations which are be solved simultaneously for S and Sc- 

45^22 7 ? 2 ^Gi 2 

2 dS ^ 2 dS 
42dE22 7?2^Gi2 

2 dSc 2 dSc 

The above equations indicate that the microdamage and transverse cracking states depend only on the strain state, 
elastic properties, and damage functions. Eurthermore, Equations ( 17 ) and ( 18 ) preclude the use of any failure criteria 
to determine the onset of matrix damage or failure. The damage functions, and slopes of the damage functions, dictate 
the manner in which the IS Vs evolve. Thus, for low strain levels, the damage IS Vs remain small and have a negligible 
effect on the response of the material. At intermediate strain levels the influence of the microdamage ISV becomes 
more prevalent. At larger strain levels, the transverse cracking ISV begins to increase more rapidly. At even higher 
strain levels, the IS Vs are large enough that the transverse and shear stiffness begin to approach zero. 


( 17 ) 

( 18 ) 
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Experimentally it has been determined that S and Sc behave as e^, thus it is convenient to introduce reduced 
damage variables Sr and Scr-^ 


Sr = S^/^ 


(19) 

( 20 ) 


'cr — >^c 




Furthermore, the use of reduced IS Vs yields polynomial damage functions. With the reduced IS Vs, the evolution 
equations. Equations (17) and (18), become 



( 21 ) 


( 22 ) 


Once Sr and Scr are determined from Equations (21) and (22), while ensuring Equations ( 8 ) and (9), the transverse 
and shear moduli can be degraded accordingly using Equations (13) and (14). The damage functions used in Equations 
(13) and (14), and the procedures for obtaining, are given in Section III. 

It is important to note that ST can account for compressive failure of a laminate resulting from fiber microbuckling 
and kink band formation within the laminae. As a material point is loaded, the fibers begin to rotate. All calculations 
are performed in the instantaneous material frame by tracking the fiber angle and transforming the field variables 
accordingly. These rotations induce large, localized shear strains which, in turn degrade the matrix shear modulus, 
allowing the fibers to rotate more easily. If the loading is compressive, this feedback will lead to a runaway instability 
and regions of microbuckled, or kinked, fibers will develop. In previous studies, where the loading was compression 
dominated, the instantaneous fiber rotation is captured in conjunction with ST. This has the advantage of predicting 
fiber kinking failure, from a knowledge of basic constituent level properties, avoiding the use of an explicit fiber 
direction compressive strength criterion, and the need for an additional coupon level test to measure compression 
strength, 

II. B. 3. Decoupling the evolution equations 

To obtain a solution for Equations (21) and (22) the damage functions must be explicitly defined. Typically fourth- 
order, or higher, polynomials are used for the microdamage and transverse cracking damage functions. The use of 
these polynomial damage functions results in two coupled, bivariate polynomial evolution equations. The solution of 
these equations produces numerous values for Sr and Scr- Only one set of solutions represents the physical behavior 
of the material, the remaining solutions are mathematical artifacts. However to determine the correct solutions, all 
solutions to the evolution equations must be obtained. This solution procedure can quickly become cumbersome and 
computationally time consuming. 

The solution of the evolution equations can be greatly simplified if the effects of microdamage and transverse 
cracking on the transverse and shear stiffness are decoupled. Microdamage primarily manifests as shear microcracks 
between fibers.^ A photograph of a [+45/ — 45] 5 composite specimen displaying these shear microcracks, taken under 
a scanning electron microscope (SEM), is presented in Figure 2. Additionally, the influence of transverse cracking on 
the shear modulus of a composite laminate is minimal.^ Thus, the evolution equations can be decoupled by assuming 
that the shear modulus is affected only by microdamage, and the transverse stiffness is affected only by transverse 


cracking. 


E 22 = E220e-c{Scr) 
G\2 = Gi20gs{Sr) 


(23) 

(24) 


Using these relationships, the decoupled evolution equations are 



7?2^ d9a{Sr) 



(25) 



(26) 


Even though the interaction between microdamage and transverse cracking is lost by decoupling the evolution equa- 
tions, the solution procedure is greatly simplified making it more amenable for implementation into a finite element 
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model where the evolution equations will need to be solved for every applicable material point at every time incre- 
ment. As stated earlier, there is evidence to suggest that the uncoupled equations for damage evolution in laminated 
composites is valid.^ 



Figure 2. SEM photograph showing shear microcracks in damaged [+45/ — 45] s specimen. 


III. Experimental Program at NASA Langley Research Center 

The response of IM7/8552 center-notched panels reinforced with buffer strip stiffeners was investigated at the 
NASA Langley Research Center (LaRC) under the supersonics program.^^ Each panel (Figure 3) consisted of a flat 
12” X 34” composite laminate. The laminate was reinforced with two 1” wide buffer strips located 1.5” from each 
edge. These buffer strips contained eight additional 0° plies not present in the rest of the laminate. Figure 4 displays 
the cross section of the stiffened regions. The stacking sequences and thicknesses for the base laminate and buffer 
strip regions are given in Table 1 . The flat plate and buffer strips were co-cured as a single part. After the pristine 
panel was cured, a centered notch that was 3” long and had notch tip diameter of 0.09735” was machined in the 
panel. 5” x 12” steel grips were adhered to the ends panel. The grips were manufactured so they would conform to 
the faces of the panel. One component of each set of grips was machined to accommodate the buffer strips, and the 
other component was left flat. The bottom grip was clamped, preventing any displacements or rotations. The top grip 
was secured to prevent any rotations or displacements in the x- or 2 ;-directions, while a tensile displacement of 5E-5 
in./sec. was applied in the ^-direction until the specimen failed catastrophically. Several strain gages, shown in Figure 
3, were affixed to the specimen and were monitored during loading of the specimen along with the applied load and 
edge displacement. 


Region 

Stacking 

Thickness (in.) 

Flat Panel 

[45°/90°/-45°/0°]5 

0.044 

Buffer Strip 

[45°/90°/-45°/0g]5 

0.088 


Table 1. Panel stacking sequence and thickness in flat and buffer strip-reinforced regions. 


The buffer strips were added to the base laminate to increase the ultimate strength of the center-notched panels 
without a substantial weight increase. Moreover, the intention of the buffer strips was to alter the failure mode of the 
panel. Additional buffer strip configurations, not presented here, were also tested. These panels exhibited self similar 
crack growth propagating from the notch tip towards the free edge. Upon reaching the buffer strips, however, the 
direction of damage propagation turns so that the damage runs parallel to the buffer strip towards the gripped edge. 
This behavior was not observed with the buffer strip configuration presented in this paper. Instead, the crack like 
damage was arrested briefly at the buffer strip boundary before proceeding through the buffer strip towards the free 
edge. Three specimens with the lay-up presented in Table 1 were tested in all. 

The elastic properties for IM7/8552 are taken from Ref. 32 and given in Table 2. The shear microdamage curve was 
obtained from a [+45/-45]5 tensile dogbone specimens. The normal transverse cracking damage curve was obtained 
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Figure 3. Schematic of center-notched, buffer strip reinforced, tensile specimen. 



Figure 4. Schematic of cross-section of buffer strip reinforced regions. 
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from multi-angle AS4/3502 specimens, presented in Ref. 2. The coefficients of this curve were scaled by the ratio of 
the elastic E 220 for IM7/8552 to the elastic E 220 for AS4/8502. 


IM7/8552 

^ci 


rIM7/8552 

^220 

rAS4/3502 

^220 


AS4/3502 
^ci 5 


i 7^0 


(27) 


This scaling procedure was used in the successful failure prediction of notched compression specimens as reported 
in Ref. 63. Figures 5 and 6 show the shear microdamage curve and normal transverse cracking damage curve as 
functions of the respective reduced IS Vs. The coefficients of the polynomial curves are presented in Table 3. 


Property 

Value 

Ell (msi) 

24.86 

E 22 (msi) 

1.78 

Z^12 

0.32 

Gi 2 (msi) 

0.68 


Table 2. Elastic properties and axial strain allowable for IM7/8552. 



Figure 5. Shear microdamage function. 



Figure 6. Normal transverse cracking damage function. 
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Function 

Co 

Cl 

C2 

C3 

C4 

C5 

9s i^r) 

1.0000 

-6.3955E-2 

-4.7716E-2 

1.0566E-2 

-8.3609E-4 

2.2997E-5 

( ^cr ) 

1.0000 

4.9809E-2 

-6.7644E-2 

-1.5854E-2 

2.2461E-3 

0.0000 


Table 3. Polynomial coefficients of damage curves. 


IV. Finite Element Model 

IV. A. Finite element discretization and boundary conditions 

One quarter of the center-notched, buffer-strip panel, presented in Section III, was discretized using ABAQUS, S4R, 
4-noded, quadrilateral and S3R, 3-noded, triangular, two-dimensional (2D), plane stress, shell elements. This mesh 
contained 9257 S4R and 42 S3R elements, with a high density of elements near the notch tips to the free edge of the 
panel. The high density regions were chosen to accurately capture the fields near the notch tip and along the expected 
axial failure path. This mesh is displayed in Figure 7. The portions highlighted in green indicate buffer strip-reinforced 
regions. 



Figure 7. Finite element mesh of quarter of buffer strip panel domain. Buffer strip reinforced elements highlighted in red. 


Displacements in the x-, y-, and 2 ;-directions as well as all rotations along the bottom edge of all the panel are 
constrained. Similarly for the top edges, the displacements in the x- and ^-directions and all rotations are prohibited. 
A uniform displacement in the ^-direction is applied dynamically at a rate of 0.886 in./sec. up to 0.0886 in. over 
a duration of 0.1 sec. Initially, mass scaling is introduced into the model on an element by element basis, using the 
ABAQUS/Explicit command ^MASS SCALING, to increase the minimum stable time increment of the problem to 
l.OE-5 sec. The input deck card ^VARIABLE MASS SCALING is also used to scale the mass of the elements every 
1000 time increments if necessary. 

IV.B. In-plane damage and failure 

Dual ISV ST is implemented into the FEM model using the commercially available ABAQUS 6.8-1 finite element soft- 
ware. The ST constitutive laws are integrated into ABAQUS/Explicit through a user material subroutine (VUMAT). 
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Axial failure is also included in the VUMAT to capture the effects of fiber breakage and other axial failure mechanisms. 
Prior to calculating the damage state of the material point, an axial maximum strain failure criterion is evaluated, 

e/ = ^ (28) 

where en is the strain in the fiber direction of the layer and is the ultimate strain. If e/ > 0 the element is deleted 
and removed from subsequent calculations. If axial failure has not occurred; then, at each integration point in an 
element obeying ST, the local strain fields are used in Equations (25) and (26) to calculate Sr and 5'cr- Equations (8) 
and (9) are enforced manually; if an ISV decreases, the value of the ISV at the end of the previous increment is used. 
Once the damage state is determined, the transverse and shear stiffness are updated using Equations (23) and (24). 
The stress field is then calculated using the degraded moduli and Equation (10). 

The IM7/8552 lamina elastic properties are given in Table 2, and the damage function coefficients are presented 
in Table 3. It should be mentioned that the axial failure criterion. Equation (28) and the degradation scheme used are 
very mesh dependent. The experimental, average, maximum value of strain near the notch (Sg-8) was used to estimate 
the maximum strain allowable for the simulation. This strain allowable was 2992 microstrain (/ie). Due to the very 
steep strain gradients near the notch, the average strain in a region the size of the element used at the notch tip in the 
simulation would be far less than the strain reported by the strain gage. Therefore, the experimentally obtained strain 
allowable was reduced by 25% to arrive at an average allowable for the axial failure strain of the elements used in 
the model: = 2250/ie. Ongoing studies will include a sensitivity analysis of to examine how it infiuences the 

overall results of the model. Additionally, future models will incorporate more sophisticated axial failure models that 
can account for pullout and friction during axial failure of a fiber composite,^^ as well as possible non-local effects,^^ 
thus, alleviating the dependence of axial failure on mesh size and density. 


V. Results 

V.A. Quantitative results compared to experimental data 

Three center-notched, buffer strip-reinforced panels were tested under uniaxial tension at NASA LaRC. Load versus 
edge displacement as well as local strain gage data was recorded. The location of the strain gages is displayed in 
Eigure 3. Each strain gage shown in the diagram represents multiple strain gages placed in equivalent locations on the 
panel. The data presented here is the average value of those analogous strain gage readings. 

A comparison of the global load-deflection curve obtained from the model to those from experiment is plotted in 
Eigure 8. The simulation is able to capture both the linear and non-linear load versus displacement response of the 
panel. There was significant variability between experimental results for the nominally same panel. The ultimate loads 
for the three panels were 37,407 lb., 30,451 lb., and 23,471 lb. Thus, the average ultimate load for the three 
tests was 30,443 lb. with a standard deviation of 6,968 lb. The finite element model yielded an ultimate load of 32,844 
which is approximately 8% greater than the average. This indicates that the 25% reduction imposed on the average 
maximum strain near the notch tip to obtain the axial failure strain allowable was slightly non-conservative, or that 
energy was dissipated via other damage mechanisms not included in the finite element model, such as delamination or 
fiber bridging. Discrepancies between the predicted ultimate load and the experimental values could also be a result 
of variability in the experiment, or a artifact of the degradation scheme used subsequent to satisfaction of the failure 
criterion. The ultimate loads are summarized in Table 4. In addition to an adequate ultimate load prediction, the 
simulation also exhibits several discontinuous load-drops, that are seen in the experimental results, prior to reaching 
the ultimate load. This may indicate that the failure and damage mechanism evolution predicted with the finite element 
model are accurate. An overview of these damage and failure mechanisms are presented in the next section. 

Eigure 9 shows the ratio of the kinetic energy (KE) to the total energy (ETotai) of the system as a function of time. 
The kinetic energy remains less than 5% of the total energy throughout the simulation except at the very beginning of 
the model (well within the linear regime), and just prior to ultimate failure. This indicates that the mass scaling and 
rapid loading rate used to speed up the total simulation time are not infiuencing the results in a non-physical manner. 

Local strain gage data for the nine gages (see Eigure 3) are plotted in Eigure 10. The linear stiffness of the numerical 
model is not in correspondence with any of the strain gage stiffness’. It should be noted that there is a discrepancy 
between the stiffness of the experimental load-displacement curve and the load-axial strain curve of Sg-1 far away 
from the notch. The reason for this discrepancy has yet to be uncovered; however, the trends predicted by the finite 
element model can still be compared to those observed in the experiment. Since the ultimate load obtained from the 
model was closest to that in Test2, local model data will be compared most closely to that test. 
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Figure 9. Kinetic energy history. 
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Load (lb.) Load (lb.) Load (lb.) 




Strain (\i£) Strain (|ie) 



Strain (|ie) 


(a) Sg-1 


(b) Sg-2 


(c) Sg-3 





(d) Sg-4 


(e) Sg-5 


(f) Sg-6 





(g) Sg-7 


(h) Sg-8 


(i) Sg-9 


Figure 10. Experimental and numerical load vs. axial strain data for buffer strip panel. 
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At Sg-1, the initial load at which the strain begins to relax in the numerical simulation is comparable to the load in 
Test2. This strain relaxation far from the notch is a result of damage occuring near the notch tip. The strain response 
at Sg-2, which is located on the buffer strip far from the notch, is linear in both the simulation and the experiment. 
Sg-3, which is located at the boundary of the buffer strip but 6 in. above the notch, exhibits an initial strain relaxation 
in both Test2 and the finite element model at the same load level as observed in Sg-1. Jumps in the strain levels at Sg-4 
and Sg-6 occur in the model at global loads significantly less than in Test2. This indicates that the damage and failure 
evolution progresses through the buffer strip more rapidly in the numerical model than in the experiment. However, 
data for Sg-6, which lies on the buffer strip directly in front of the notch, matches well with the experiment. The 
same can be seen at Sg-7. Experimental data at the notch tip (Sg-8) displays a high level of non-linearity. This is 
not observed in the FE model and may be a facet of the extremely steep strain gradients at the notch tip which would 
require an extremely fine mesh at that region. Similarly, the numerical model does not capture the non-linearity shown 
in the experiment at Sg-9 in the center of the panel directly above the notch. A relatively large element is present in 
the model at that location and is used to represent that strain gage; a more refined mesh in that region may capture the 
appropriate response, and this is currently being studied. Additionally, poor performance of the FEM model at Sg-8 
and -9 may be due to the lack of transverse- shear damage coupling. 

V.B. Analysis damage and failure mechanisms observed in numerical simulation 

The results of the finite element model indicate that large degrees of microdamage, transverse cracking, and axial 
failure are confined to regions in front of the notch tip. Thus, Sg-1, Sg-2, and Sg-3 (see Figure 3), which lie in regions 
far from the notch, should not be directly influenced by damage or failure, but rather are affected through redistribution 
of the strain fields. Figure 11 shows the evolution of the microdamage and transverse crack IS Vs with respect to the 
applied displacement, in each of the top four plies, at the strain gage locations near and in front of the notch tip (Sg- 
4-Sg-9). The ISVs are normalized by their maximum values = 12.08 psi^/^ and 5'^"^ = 3.52 psi^/^) so that 

a value of one in Figure 1 1 corresponds to complete degradation of the appropriate modulus (see Equations (23) and 
(24)). 

Similar damage evolution can be seen at Sg-4, Sg-5, and Sg-6. Microdamage evolves progressively, and equally, 
in the -\-45° and -45^ plies. Very low levels of Sr arise at these gages in the 0° plies, and no microdamage is observed 
in the 90° ply until catastrophic failure. Transverse crack evolution remains dormant in all plies, until after a 0.13 in. 
edge displacement is achieved; after which the degree of damage due to transverse cracking rises rapidly in the 90° 
ply until catastrophic failure. This increase in Scr is associated with propagation of axial failure from the notch tip, 
as is evidenced in the proceeding damage contour analysis. In the -f 45°, -45°, and 0° plies, Scr remains zero until 
catastrophic failure. 

At the strain gage on the boundary of the buffer strip closest to the notch, only very low levels of microdamage 
are observed in the +45° and -45° plies. However, both Sr and Scr jump suddenly at an edge displacement of 0.1 13”. 
This jump in the ISVs is associated with axial failure evolving from the notch meeting the buffer strip boundary. 
At the notch tip (Sg-8), transverse cracking in the +45° and 90° layers develop rapidly at low displacement levels. 
Expeditious advancement in Scr is accompanied by a steady increase of Sr in the +45° and -45° layers near the notch 
tip. Less severe microdamage evolution also occurs in the 90° and 0° laminae. At a global displacement of 0.1”. there 
is a large jump in Scr in the -45° layer. Scr remains zero in the 0° plies; this is most likely a result of axial failure of 
these layers at the notch tip before transverse cracks can develop. Above the notch at Sg-9, low levels of microdamage 
in the +45° and -45° are observed. Overall, the evolution of microdamage in the plies is progressive; whereas, the 
advancement of transverse cracks is sudden and discontinuous. 

Figures 12-15 show matrix microdamage and transverse cracking contours, as well as axial failure paths, at several 



Load (lb.) 

Deviation from 
Experimental Average 

Testl 

37,407 

23% 

Test2 

30,451 

0.0003% 

Test3 

23,471 

23% 

Average 

30,443 

N/A 

FEM 

32,884 

8% 


Table 4. Ultimate load information for experiments and finite element model. 
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Normalized Damage Variable S Normalized Damage Variable S 



(a) Sg-4 





(b) Sg-5 




(c) Sg-6 



(d) Sg-7 


(e) Sg-8 


(f) Sg-9 


Figure 11. Microdamage and transverse crack evolution at strain gage locations in front notch tip. 
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(g) Sr, t = 0.0054667s 


(h) Sr, t = 0.0061334s 


(i) Sr, t = 0.0064001s 



(j)S'cr,l = 0.0054667s (k) Scr, * = 0.0061334s (1) Scr, * = 0.0064001s 

Figure 12. Microdamage, transverse crack evolution and axial failure in top 45^ ply. 
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(j) Scr,t = 0.0054667s 


(k) S'cr, t = 0.0061334s 


(1) Scr, t = 0.0064001s 


Figure 13. Microdamage, transverse crack evolution and axial failure in top 90° ply. 
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(g) Sr, t = 0.0054667s (h) Sr, t = 0.0061334s (i) Sr, t = 0.0064001s 



(j) Scr, t = 0.0054667s (k) Scr, t = 0.0061334s (1) Scr, t = 0.0064001s 

Figure 14. Microdamage, transverse crack evolution and axial failure in top -45° ply. 
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(j) Scr, t = 0.0054667s (k) Scr, t = 0.0061334s (1) Scr, t = 0.0064001s 

Figure 15. Microdamage, transverse crack evolution and axial failure in 0° plies. 
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key time points in the top four layers (45°, -45°, 90°, and 0°) of the panel acreage and buffer strip. These snapshots 
depict the evolution of, as well as the interactions between, the three damage mechanisms modeled in this work. The 
damage variables, Sr and Scr are represented continuously in the colored contours. Blue indicates zero, or very little, 
damage, while red indicates the maximum damage (i.e. when the stiffness has degraded to zero). Turquoise, green, 
yellow, and orange bridge the scale between minimum damage (blue) and maximum damage (red). Elements that have 
satisfied the axial failure criterion, or exceed some excessively large strain value, are considered to not contribute any 
stiffness to the panel and are removed. Elements in the buffer strip region are outlined in green to display the effects 
of the pad-up in this region on the crack path and damage evolution. 

In the outer 45° layer, prior to initiation of axial failure, a very small amount of microdamage is present near the 
notch tip and is accompanied by a localized region of maximum transverse cracking damage. Shortly after this, axial 
failure occurs in these elements. When the maximum strain criterion is satisfied in an element, that element is deleted, 
and the corresponding void that develops essentially represents a macroscopic crack in the material. Additionally, any 
elements that have reached and have a transverse strain that has surpassed the axial strain allowable, are also 

removed to prevent excessive element distortion in the simulation. This is a reasonable assumption since transverse 
failure strains are typically lower than axial failure strains, and an element that has reached has most likely been 
loaded predominantly in the transverse direction and is not contributing axially. When the crack reaches the buffer 
strips it changes direction and runs parallel to the buffer strip. The crack is then arrested for a brief moment before 
it advances laterally into the buffer strip (Eigure 12(c)). Once inside the buffer strip, the crack turns again and runs 
along the buffer strip. The subsequent point of arrest for the crack is in a transition region between a very dense mesh 
and a coarser mesh. It may be that this transition region is influencing the behavior of the crack in a non-physical 
way; however, similar behavior was observed in regions with a uniformly dense mesh (see Eigures 12(c) and 12(f)). 
A crack branches from the original crack and propagates perpendicular to the length of the buffer strip. Shortly after 
this branched crack forms, numerous other branched cracks advance through the buffer strip leading to catastrophic 
failure of the panel. The evolution of the crack in the 45° layer is not accompanied by much microdamage. However, 
large regions exhibiting high levels of transverse cracking surround the macroscopic crack is evident. 

The same macroscopic crack progression is observed in the top 90°, -45°, and 0° layers. Additionally, the accu- 
mulation in all four plies is very similar. The 0° and 90° plies exhibit more microdamage encompassing a larger area 
surrounding the crack as it advances. The 90° layer exhibits widespread damage due to transverse cracking (Eigure 
13. In the 0° layer, high values of Scr are limited to the vicinity surrounding the axial failure path. 

The individual damage mechanisms (i.e. microdamage, transverse cracking, and axial failure) are seen to be 
interactive. Accumulation of Sr and Scr near the notch tip leads to a local transfer of load into the 0° ply from 
adjacent layers, which subsequently leads to local axial failure. Similarly, the fields surrounding the advancing crack 
lead to rapid matrix degradation and further crack advancement. Upon, reaching the buffer strip, the path of least 
resistance to axial failure lies along the longitudinal direction of the buffer strip. However, substantial regions of 
severe matrix degradation in the non-0° plies in the buffer strip provide an avenue for axial failure to occur within 
the 0° layers of the buffer strip through load transfer. Once the reinforcement provided by the buffer strip begins to 
diminish, catastrophic failure is imminent and this is captured in the present work. 

VI. Conclusions 

A finite element based implementation of progressive damage accumulation and failure, or in short, progressive 
failure analysis (PEA), that is capable of capturing three predominant damage mechanisms in an ERL (namely: matrix 
shear microdamage, matrix transverse cracking, and axial failure) has been presented. The matrix damage mechanisms 
are modeled using a dual ISV formulation of ST, which accounts for the damage development mechanisms via the 
energy required to create structural changes associated with the damage modes in the material. A maximum strain 
cut-off is chosen to capture axial failure within layers of a laminate. The capabilities of the model are demonstrated 
by analyzing a center-notched, buffer-strip reinforced composite panel, for which there is experimental data. Good 
quantitative agreement with global load-displacement, and local strain gage data is exhibited. Moreover, the result- 
ing damage evolution in the model provides insight into the interactions between the individual damage mechanisms 
within the composite that lead to the global observed failure modes. To obtain subtle details, transparency, and acces- 
sibility into the evolution of damage modes in a complex composite structure with experimentation alone would be 
difficult. 

While a significant number of important damage accumulation and failure mechanisms have been captured, de- 
lamination failure has not been included in the model, but it can be incorporated by including DCZM elements at the 
interfaces between adjacent shell or three-dimensional sublaminate layers. Eurther investigations will focus on the 
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influence of interlamina damage on the response of the buffer strip panel. Additionally, the axial failure theory used 
does not properly capture the physics involved in axial damage of FRLs. A sensitivity analysis will be performed to 
determine the influence of the axial failure strain on the ultimate load of, and failure path within, the panel. Future 
efforts will also attempt to forgo the maximum strain criterion for axial failure in favor of a physically-based, non-local 
axial damage theory,^^ or a complete micromechanics simulation of axial failure that includes pull-out and frictional 
slip.^^ 
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